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1 Introduction and Notation 

These notes attempt to present the essential ideas of configuration interaction 
(CI) theory in a fairly detailed mathematical framework. Of all the ab initio 
methods, CI is probably the easiest to understand — and perhaps one of the hardest 
to implement efficiently on a computer! The next two sections explain what 
the CI method is: the matrix formulation of the Schrodinger equation H^f = 
E^f . The remaining sections describe various simplifications, approximations, 
and computational techniquies. 

I have attempted to use a uniform notation throughout these notes. Much 
of the notation is consistent with that of Szabo and Ostlund, Modern Quantum 
Chemistry [1]. Below are listed several of the commonly- used symbols and their 
meanings. 

TV The number of electrons in the system. 

n a The number of alpha electrons. 

rip The number of beta electrons. 

n The number of orbitals in the one-particle basis set. 

5ij Kronecker delta function, equal to one if % = j and zero otherwise. 

H The exact nonrelativistic Hamiltonian operator. 

H The Hamiltonian matrix, i.e. the matrix form of H, in whatever TV-electron 
basis is currently being used. 

Hij The i, j-th element of H, equal to ($>i\H\$j), where ^ and $j are TV-electron 
CI basis functions. 

Xj The space and spin coordinates of particle i. 

Yi The spatial coordinates of particle %. 
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4>i The z-th one-particle basis function (orbital). Usually denotes a spin-orbital 
obtained from a Hartree-Fock procedure. May also be written simply as z. 

Xi The z-th one-particle basis function (orbital). Usually denotes an atomic spin- 
orbital. 

The z-th A^-electron basis function. Usually denotes a single Slater determi- 
nant, but may also be a configuration state function (CSF). 

Usually denotes an eigenfunction of H. The exact nonrelativistic wavefunc- 
tion if a complete basis is used in the expansion of H. 

\<l> r a ) An iV-electron basis function which differs from some reference function |$o) 
by the replacement of spin-orbital a by spin-orbital r. Usually implies a single 
Slater determinant. 

\ab . . . c) A Slater determinant with spin-orbitals a, b, . . . c occupied, i.e. 

a (xi) 0fc(xi) ... c (xi) 
a (x 2 ) 6 (X 2 ) • • • C (X 2 ) 

a (x7v) ^(x^v) ... c (Xiv) 

(i\h\j) One-electron integral in physicists' notation (i and j are spin-orbitals). 
More explicitly, this is 

(i\h\j) = J 0*(x 1 )/ i (x 1 )0 7 (x 1 )dx 1 (1.2) 

[i\h\j] One-electron integral in chemists' notation, where z and j are spin-orbitals. 
Equivalent to (i\h\j). 

(i\h\j) One-electron integral in chemists' notation (i and j are spatial orbitals). 

(ij\\kl) Antisymmetrized two-electron integral, equal to (ij\kl) — (ij\lk). 



\(j) a (j) b . . . (j) c ) = 
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(ij\kl) A simple two-electron integral, in physicists' notation, where z, j, k, and / 
are spin-orbitals. This is 

r 1 

(ij\kl) = / 0*(xi)0*(x 2 ) — fc (xi)0/(x 2 )dxidx2 (1.3) 

[ zj|&;/ ] A simple two-electron integral in chemists' notation, where z, j, k, and 
/ are spin-orbitals. This is 

[ij\kl] = \ ^(xi)^-(xi)— 0J(x 2 )0z(x 2 )dxidx 2 (1-4) 

J 1*12 

(r/|/c/) A simple two-electron in chemists' notation where z, j, /c, and / are spatial 
orbitals. This is 

fe'lfeO = / #(n)0 7 -(n)— ^(r2)0/(r2)^^r 2 (1.5) 

J 1*12 

aj Second-quantized creation operator for orbital z. 
Second-quantized annihilation operator for orbital i. 
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2 Fundamental Concepts 
2.1 Scope of the Method 

Configuration interaction (CI) is a method for solving the nonrelativistic Schrodinger 
equation 

£*(r,R) = (5;^v^ + 5;k 2 + £ l^+E^ + E^tMl^^R) 

I A i 1 A>B K AB Ai r Ai i>j Tij J 

(2.6) 

where i,j denote electrons and A, B denote nuclei, with = — rj\, Rai = 
|Ra— r*i | , and Rab = \Ra~ Rb|- Typical applications of the CI method employ the 
Born-Oppenheimer approximation, whereby the the motions of the electrons are 
treated as uncoupled from those of the nuclei. Thus the "electronic" Shrodinger 
equation is solved at discrete sets of fixed nuclear positions 



£ e *e(r;R) = 




The Born-Oppenheimer approximation is invoked so often in computational quan- 
tum chemistry that the subscripts in the preceeding equation are usually sup- 
pressed and the equation is written simply as = E^>. However, it is important 
to remember that the electronic energy E e is an artifact of the Born-Oppenheimer 
approximation and is not as physically meaningful as the total energy of a system. 
Within the Born-Oppenheimer approximation, we estimate the total energy by 
adding the nuclear-nuclear repulsion energy and the nuclear kinetic energy to the 
total electronic energy E e of equation (2.7). 

While the CI method can be extended to incorporate some relativistic effects 
(e.g. spin-orbit terms), this is not generally done; these notes will be concerned 
only with the nonrelativistic Hamiltonian (2.7). 
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2.2 Why Configuration Interaction? 

In the first paper on quantum mechanics, Heisenberg used matrix mechanics to cal- 
culate the frequencies and intensities of spectral lines [2]. Later, when Schrodinger 
discovered wave mechanics, it was quickly shown that the Schrodinger and Heisen- 
berg approaches are mathematically equivalent [3, 4]. Given the ease with which 
matrices may be implemented on a computer, it is entirely natural to attempt 
to solve the molecular time-independent Schrodinger equation H^> = E^> using 
matrix mechanics. 

Matrix mechanics requires that we choose a vector space for the expansion of 
the problem. For the case of an TV-electron molecule, our wavefunction must be 
expanded in a basis of iV-particle functions (the nuclei need not be considered in 
the electronic wavefunction, if we have invoked the Born-Oppenheimer approx- 
imation). How do we construct the iV-particle basis functions? Here we follow 
the arguments of Szabo and Ostlund [1], p. 60. Assume we have a complete set 
of functions {x^i)} °f a single variable X\. Then any arbitrary function of that 
variable can be expanded exactly as 

$Oi) = y Ea i Xi{xi). (2.1) 

i 

How can we expand a function of two variables X\ and x 2 which have the same 
domain? If we hold x 2 fixed, then 

$(x u x 2 ) = Y, a ii x 2)Xi{xi)- ( 2 - 2 ) 

i 

Now note that each expansion coefficient di(x 2 ) is a function of a single variable, 
which can be expanded as 

a i( x 2) = YsKXjfa)- (2.3) 

j 

Substituting this expression into the one for <fr(xi,x 2 ), we now have 



(2.4) 
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a process which can obviously be extended for x 2 , • • • , %n)- 

Let us now collect the spin and space coordinates of an electron into a variable 
x. We can write a spin orbital as x( x )- The result analogous to equation (2.4) 
for a system of N electrons is 

$(xi,x 2 , . . . ,xjv) = E &y\..wXi(xi)x/(x 2 ) . . .Xn(*n) (2.5) 

ij...N 

However, the wavefunction must be antisymmetric with respect to the exchange 
of the coordinates of any two electrons 1 For the two-particle case, the requirement 

$(xi,x 2 ) = -$(x 2 ,xi) (2.6) 
implies that 6^ = — bji and ba = 0, or 

$(xi,x 2 ) = E & u'[Xi( x i)Xi( x 2) -Xj(xi)Xi(x 2 )] 

= E2 1/2 6,,|x^} (2.7) 

More generally, an arbitrary iV-electron wavefunction can be expressed exactly as 
a linear combination of all possible iV-electron Slater determinants formed from 
a complete set of spin orbitals {xi(x)}- If we solve the matrix mechanics problem 
H|\I/} = E\fy) in a complete basis of iV-electron functions as just described, we 
will obtain all electronic eigenstates of the system exactly. If our iV-electron basis 
functions are denoted the eigenvectors of H are given as 

|tf;>=Eci;l*<> (2.8) 

i 

if there are / possible A^-electron basis functions (I will be infinite if we actually 
have a complete set of one electron functions Xi)- The matrix H is constructed 
so that Hij = (<$>i\H\<&j) for i,j = 1,2, ... ,1. The matrix elements Hij may be 
written in terms of one- and two-electron integrals according to "Slater's rules," 
as discussed in section 2.4. 

1 According to the antisymmetry principle for fcrmions, of which the Pauli principle is a direct result. 



2 FUNDAMENTAL CONCEPTS 



10 



The iV-electron basis functions can be written as substitutions or "excita- 
tions" from the Hartree-Fock "reference" determinant, i.e. 

|*)=c |$o) + EcK)+ E cZ\<f>2)+ E &^ c ) + ... (2.9) 

ra a<b,r<s r<s<t,a<b<c 

where \§ r a ) means the Slater determinant formed by replacing spin-orbital a in |$o) 
with spin orbital r, etc. Every TV-electron Slater determinant can be described 
by the set of N spin orbitals from which it is formed, and this set of orbital 
occupancies is often referred to as a "configuration." Thus the "configuration 
interaction" method is, in its most straigtforward implementation, nothing more 
or less than the matrix mechanics solution of the time-independent non-relativistic 
electronic Schrodinger equation H^> = E^f. One of the great strengths of the 
CI method is its generality; the formalism applies to excited states, to open- 
shell systems, and to systems far from their equilibrium geometries. By contrast, 
traditional single-reference perturbation theory and coupled-cluster approaches 
generally assume that the reference configuration is dominant, and they may fail 
when it is not. 

In practice, one does not have a complete set of one-particle basis functions 
{Xi(x)}; typically one assumes that the incomplete one-electron basis set is large 
enough to give useful results, and the CI procedure is not modified. The quality of 
the one-particle basis set can be checked by comparing the results of calculations 
using progressively larger basis sets. 

It is also possible to reduce the size of the iV-electron basis set. If we desire 
only wavefunctions of a given spin and/or spatial symmetry, as is usually the case, 
we need include only those iV-electron basis functions of that symmetry, since the 
Hamiltonian matrix is block-diagonal according to space and spin symmetries. 
This point is discussed further in section 4.1. If one performs the matrix mechanics 
calculation using a given set of one-particle functions {xi(x)} and all possible N- 
electron basis functions (possibly symmetry-restricted), the procedure is 

called "full CI." The full CI corresponds to solving Schrodinger's equation exactly 
within the space spanned by the specified one-electron basis. If the one-electron 
basis is complete (it never is in practice, but it may be in theory), then the 
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procedure is called a "complete CI" [5]. 

Unfortunately, even with an incomplete one-electron basis, a full CI is compu- 
tationally intractable for any but the smallest systems, due to the vast number 
of TV-electron basis functions required (the size of the CI space is discussed in 
section 4.4). The CI space must be reduced somehow — hopefully in such a way 
that the approximate CI wavefunction and energy are as close as possible to the 
exact values. The effective reduction of the CI space is a major concern in CI 
theory, and we will discuss some of the more popular approaches in these notes. 

By far the most common CI approximation is the truncation of the CI space 
expansion according to excitation level relative to the reference state (equation 
2.9). The widely-employed CI singles and doubles (CISD) wavefunction includes 
only those iV-electron basis functions which represent single or double excitations 
relative to the reference state. Since the Hamiltonian operator includes only one- 
and two-electron terms, only singly and doubly excited configurations can interact 
directly with the reference, and they typically account for about 95% of the corre- 
lation energy in small molecules at their equilibrium geometries [6] . Truncation of 
the CI space according to excitation class is discussed more thoroughly in section 
4.2. 

2.3 The Correlation Energy 

Approximate CI methods can be judged according to what fraction of the cor- 
relation energy they recover. The correlation energy is defined as the difference 
between the energy in the Hartree-Fock limit (Ejjp) and the exact nonrelativistic 
energy of a system (So) 

E covv = Sq — Ehf (2.10) 

This energy will always be negative because the Hartree-Fock energy is an upper 
bound to the exact energy (this is guaranteed by the variational theorem, as 
explained in section 3). The exact nonrelativistic energy So could, in principle, 
be calculated by performing a full CI in a complete one-electron basis set. If we 
have an incomplete one-electron basis set, then we can only compute the basis set 
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correlation energy, which is the correlation energy for a given one-electron basis. 
For convenience, the basis set correlation energy is often simply referred to as the 
correlation energy. 

The correlation energy is the energy recovered by fully allowing the electrons 
to avoid each other; Hartree-Fock improperly treats interelectron repulsions in 
an averaged way. 2 However, there is some inconsistency in this line of thinking. 
When a molecule is pulled apart, the electrons shouldn't need to avoid each other 
as much, so the magnitude of the correlation energy should decrease. In fact, the 
opposite is true, as shown by the basis set correlation energies given in Table 1 
for H2O at three different geometries. 

Table 1: Electron Correlation in H2O with a DZ Basis. 



Geometry 


E corr (hartree) a 


R e 


-0.148028 


1.5 R e 


-0.210992 


2.0 R e 


-0.310067 



a Data from reference [6] . 



The correlation energy increases at stretched geometries, because our defini- 
tion of the correlation energy in equation (2.10) includes not only the concept of 
electrons avoiding each other, which is called the "dynamical" correlation energy, 
but also a more subtle effect called the "nondynamical," or "static" correlation 
energy. Nondynamical correlation energy reflects the inadequacy of a single refer- 
ence in describing a given molecular state, and is due to nearly degenerate states 
or rearrangement of electrons within partially filled shells. Shavitt [7] has pointed 
out this deficiency in the definition of the correlation energy, and has suggested 
that perhaps a multiconfigurational Hartree-Fock method may be more useful in 
the definition of correlation energy. 

2 Some texts talk about the "average" electron repulsion term in the Fock operator; I find this misleading in 
that Hartree-Fock uses the same instantaneous interelectron repulsion term as CI — it's the same Hamiltonian! The 
restriction to a single Slater determinant is what causes the averaging of interelectron repulsions. 
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Siegbahn [8] offers the following explanation of the difference between dynam- 
ical and nondynamical correlation energies: 

"In many situations it is further convenient to subdivide the correla- 
tion energy into two parts with different physical origins. For chemical 
reactions where bonds are broken and formed, and for most excited 
states, the major part of the correlation energy can be obtained by 
adding only a few extra configurations besides the Hartree-Fock config- 
uration. This part of the correlation energy is due to near degeneracy 
between different configurations and has its origin quite often in artifacts 
of the Hartree-Fock approximation. The physical origin of the second 
part of the correlation energy is the dynamical correlation of the motion 
of the electrons and is therefore sometimes called the dynamical corre- 
lation energy. Since the Hamiltonian operator contains only one- and 
two-particle operators this part of the correlation energy can be very 
well described by single and double replacements from the leading, near 
degenerate, reference configurations." 

2.4 Slater's Rules 

Whether we perform a full CI or only a limited CI, we must be able to express H 
in matrix form so that we can diagonalize it and obtain the eigenvectors and eigen- 
values of interest. In this section we discuss Slater's rules (or the Slater-Condon 
rules [9, ?, 10]), which allow us to express matrix elements Hij = ($>i\H\$j) in 
terms of one- and two-electron integrals. At the moment, we will express these 
results in terms spin-orbitals using physicist's notation. The one-electron integrals 
are written as 




(2.11) 



and the two-electron integrals are written as 



(ij\\kl) = (ij\kl) - (ij\lk) 



(2.12) 
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where 

(ij\kl) = 1 0*( ri )0*(r 2 )— M^M^drrfri (2.13) 

Before Slater's rules can be used, the two Slater determinants must be arranged 
in maximum coincidence. Remember that switching columns in a determinant 
introduces a minus sign. For instance, to calculate (<$>i\H\<$>2), where we have 

|$i) = \abcd) (2.14) 

|$ 2 ) = \crds) (2.15) 

then we must first interchange columns of |$i) or \<&2) to make the two deter- 
minants look as much alike as possible. For example, we may rearrange ($2) 
as 

^2} = \crds) = —\crsd) = \srcd) (2.16) 

After the determinants are in maximum coincidence, we see how many spin or- 
bitals they differ by, and we then use the following rules: 

1. Identical Determinants: If the determinants are identical, then 

N N 

= J2( m \h\ m ) + Yl (mn\\mn) ( 2 -17) 

m m>n 

2. Determinants that Differ by One Spin Orbital: 

= \... mn ...) (2.18) 
I $2) = \'-'pn--) 

($i|if|$ 2 ) = (m\h\p) + Y,( mn \\P n ) 

n 

3. Determinants that Differ by Two Spin Orbitals: 

= \... mn ...) (2.19) 
|$ 2 ) = \...pq...) 
($ 1 \H\$ 2 ) = {mn\\pq) 
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4. Determinants that differ by More than Two Spin Orbitals: 



= \---mno---) (2.20) 
1^2} = \---pqr---) 



<S>i\H\<S> 2 ) = 



The derivation of these rules can be found in Szabo and Ostlund [1], section 
2.3.4 (pp. 74-81). 



3 THE VARIATIONAL THEOREM 



16 



3 The Variational Theorem 
3.1 The Method of Linear Variations 

In this section we show that the method of linear variations (also called the Ritz 
method [11]) is equivalent to the matrix formulation He = Ec of the Schrodinger 
equation H\^) = E\ty). Our treatment is similar to that of Szabo and Ostlund 
[1], p. 116. The linear variation method states that, given the linear expansion 

W=Eci|^> (3.1) 

i 

we vary the coefficients q so that we minimize E = (^\H\^) / {^f\^) . We begin by 
requiring that the wavefunction be normalized so that (^|^) = 1. Normalization 
means that we cannot minimize E simply by solving 

-^-{V\H\V) = fc = l,2,...,iV (3.2) 

because the q's are not independent. In this case we have a constrained min- 
imization, so we apply Lagrange's method of undetermined multipliers, and we 
minimize the functional 

C = (#|#|#) - E((V\V) - 1) (3.3) 

which has the same minimum as E when is normalized. When we substitute 
equation (3.1) into equation (3.3), we obtain 

£ = Ec*9<<N# ft) " E (Ec*^!^.) - l) (3.4) 

ij \ ij / 

which we may rewrite as 

£ = E dtcjHy - E (e dfaSu - l) (3.5) 
ij \ ij ) 

where of course Hij = ($>i\H\$j) and Sij = ($>i\$j). Now set the first variation in 
C equal to zero: 

SC = E SdicjHij - EJ2 SclcjSij + E ({ScjHij - E E = (3.6) 

ij ij ij ij 
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Since the summations run over all i and j, and since Hij = H*- and Sij = S*^ we 
can simplify to 



Since each term is a sum of a number and its complex conjugate, the imaginary 
parts will cancel. However, the real part will not necessarily be zero; in fact, since 
all the dci's are arbitrary (that is the whole point of using Lagrange's method), 
then for 5C to be zero, the term in brackets must be zero. We may rewrite this 
condition as a matrix equation 



If the basis functions { } are chosen orthonormal (as is usually the case), then 
S = I the identity matrix, and we have He = Ec. Of course c is the column-vector 
representation of in the basis { }. We thus have two equivalent ways of 
viewing a CI — either as the matrix formulation of the Schrodinger equation within 
the given linear vector space of TV-electron basis functions, or as the minimization 
of the energy with respect to the linear expansion coefficients q of (3.1), subject to 
the constraint that the wavefunction remain normalized. Another way of viewing 
the results of this section is to note that only eigenvectors of the Hamiltonian 
matrix H are stable with respect to variations in the linear expansion coefficients. 

At this point it is reasonable to ask why we wish to minimize the energy by 
varying the coefficients in equation (3.1). How do we know that this will give us 
the best estimate of the wavefunction? There are two answers to this. First, as 
we have just shown, minimizing the energy by variation of the linear expansion 
coefficients gives the Schrodinger equation in matrix form; thus the procedure is 
justified a posteriori by the validity of its result. The other reason is that, for the 
ground state, the linear expansion in equation (3.1) gives an expectation value for 
the energy E which is always an upper bound to the exact nonrelativistic ground 
state energy Sq. We will prove this assertion in the next section; the result is 
called the Variational Theorem. The best estimate of E, then, is the minimum 
value which can be obtained by varying the coefficients in equation (3.1) (while 



5C = $ c i Hij c j ~ ESijCj + complex conj. = 

i L 3 



(3.7) 



He = ESc 



(3.8) 
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also maintaining normalization). These arguments also hold for excited states, so 
long as each excited state is made orthogonal to all lower states. 



3.2 Variational Theorem for the Ground State 

One particularly nice feature of the CI method is that the calculated lowest en- 
ergy eigenvalue is always an upper bound to the exact ground state energy. Our 
approximate wavefunction |0) can always be expressed as a linear combination 
of the exact nonrelativistic eigenvectors which span the entire iV-electron 

space. 

|e> = Ecil*i> (3.9) 

i 

and the energy is given by 

F <e|g|e) 

Now if |0) is normalized, and we substitue the expansion over exact eigenfunctions 
(equation 3.9) into the equation above, we obtain 

£ = £c*c^ (3.11) 

i 

where Si is the zth energy eigenvalue, i.e. H\^i) = S^j). Now subtract So, the 
exact nonrelativistic ground state energy, from both sides to obtain 

E-S = E4ciSi-S (3.12) 

i 

or 

£-£o = £c*c*(£-£o) (3.13) 

i 

since |G) is normalized and EiC*Q = 1. Since Si are greater than or equal to 
Sq for all values of i and the coefficients c*Ci are necessarily non-negative, the 
right hand side of equation (3.12) is non-negative, so that E > So. We should 
also point out that the variational theorem holds for the Hartree-Fock method 
as well as for CI, since equation (3.10) is valid for the Hartree-Fock energy — for 
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a given set of MO's, the HF energy can be formulated as a (trivial) 1 x 1 CI 
eigenvalue problem. In a similar manner, the MCSCF method (where MO's and 
CI coefficients are optimized) is also "variational." 

It should be clear that instead of using the exact nonrelativistic eigenfunctions 
in equation (3.9), we could also have used an expansion over the exact eigen- 
functions within the one-electron space spanned by |0) (i.e. we could expand the 
approximate CI wavefunction |0) in terms of the full CI wavefunctions) . This 
means that not only is the approximate energy an upper bound to the exact non- 
relativistic ground-state energy, but it is also an upper bound to the full CI energy 
in the given one-electron basis. 

3.3 Why are Coupled-Cluster and MBPT Energies not Variational? 

Electron correlation methods other than CI may not be variational. For example, 
consider the coupled-cluster energy expression 



E = 



(So I 


e~ T He T \ 


S ) 


(S | 


S ) 



(3.14) 



If the operator e T is not trunctated, then we know that e r |<I>o) = l^o)- Generally, 
however, the operator is truncated. Let us define e T '|$o) = |@a) for our truncated 
T' . Now define (<I>o|e~ T ' = {®b\- Note that in general |0#) ^ |0a), which would 
have occured had we used ($o|(e T )^ cm the left. Then the energy expression is 




(3.15) 



(3.16) 



(3.17) 



(Sol^o) 
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At this point we can go no farther, because the terms c\d{ may be negative, in 
contrast to the situation in equation (3.12). 

For completeness, we also show that MBPT energies are not variational. The 
nth order MBPT wavefunction may be written [12] as 



k 

/ \ « ■■■■ ■ " ■ ■ 

l©Ltp T > = \*o) + E 

k=l 



V(l - l^oX^ol) 
Eq — H 



$o)l (3.18) 



where the sum is over "linked diagrams" only. The nth order energy is then given 

by 

E { mLt = <*o|#|efe&) (3.19) 

Since this integral is not symmetric, the energy is not variational. Only the 
first-order perturbation theory energy (which is also the Hartree-Fock energy) is 

variational, since it uses |©mbpt) = l^o)- 



3.4 Application of the Variational Theorem to Other States 

In this section, we parallel the arguments of Pauling and Wilson [13], p. 186. So 
far, we have shown only that the energy calculated as the expectation value of 
some trial function must be an upper bound to the true ground state energy So. 
In certain cases, we may derive a similar result for other states. If we take a trial 
function |$o) such that the first k coefficients in equation (3.9) are zero, then we 
may subtract Sk from equation (3.11) to obtain 

E-e k = Y J c\c i {e i -Sk) (3.20) 

i 

Now since we've assumed Cj = 0, j = 0, 1, . . . , k, this simplifies to 

E-S k = Y, ~ Sh) (3-21) 

i=k 

Once again, we can see that every term on the right side is nonnegative, so E—Sk > 
0. 
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There are any number of cases in which we have a trial function of the form just 
described. Consider, for example, a calculation on a triplet state for a molecule 
which has a singlet ground state. If our trial function is constrained to be a triplet, 
then all singlet eigenfunctions |\E^} will have zero coefficients in the expansion of 
the trial function. In this case, the energy we minimize from the triplet trial 
function will be an upper bound to the lowest triplet energy, even though there is 
a lower-lying singlet state. Similar arguments can be made for spatial symmetry. 

3.5 Convergence of the Wavef unction 

A consequence of the variational theorem is that as the energy E of an approximate 
variational wavefunction approaches the exact energy So, the approximate wave- 
function |$) approaches the exact one |^o}- This follows from equation (3.12), 
which shows that as the energy E is minimized (or, equivalently, as E — So is 
minimized), then T,iC*Ci(Si — So) is minimized; that is, the sum of squares of the 
absolute values of the coefficients of excited states with weight factors (Si — So) 
is minimized. It is apparent that these weight factors might not be optimal if we 
want the |$) which gives the best value for a property other than the energy, such 
as dipole moment. However, in the limit that E is minimized with a sufficiently 
large basis so that E = So, then T,i=i c*Ci(Si — So) = 0, or cq = 1, implying that 
|$) = |^o}- Once we have the exact wavefunction, then of course all properties 
can be computed exactly (again, within the Born-Oppernheimer approximation 
and neglecting relativistic effects). To summarize, variational improvements in 
the energy give improvements in the approximate wavefunction, which in turn 
improves the values of all other properties; however, these other properties will 
not necessarily converge as fast as the energy with respect to TV-electron basis set 
improvement. 
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3.6 Variational Theorem Bounds on Excited States 

Just as the lowest eigenvalue has been shown to be an upper bound to the exact 
ground-state energy, more generally, any eigenvalue E{ can be shown to be an 
upper bound to the corresponding exact excited state energy S{ [14]. In fact, 
one can also show that as other iV-electron basis functions are added to the CI 
procedure, the eigenvalues obey the MacDonald-Hylleraas-Undheim relations [14, 



15] 




(m+1) 



(to) 



(3.22) 



where m is the number of iV-electron basis functions. 
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4 Reducing the Size of the CI Space 



4.1 Symmetry Restrictions on the CI Space 

In this section we discuss some general considerations concerning which iV-electron 
basis functions should be included in the CI space (given that, in the general case, 
there are too many for us to include all of them). Certainly if we can find a class 
of iV-electron functions which rigorously have zero Hamiltonian matrix elements 
with the desired CI wavefunction, then none of these basis functions will contribute 
at all to our approximate wavefunction and they should not be included in the CI 
space; the Hamiltonian will be block diagonal, and all the functions in this class 
will contribute to the wrong block. We will now prove the following: consider an 
operator A which commutes with H . If 



and 
then 



=ai|$i) 
i|$ 2 ) = a 2 |$2>,ai ^ a 2 



(4.1) 
(4.2) 



($i|#|$ 2 ) =0 (4.3) 
First we show that H\& 2 ) is an eigenfunction of A with eigenvalue a 2 . Define 

H\<S> 2 ) = |$ 2 ) (4.4) 



Now apply A to |$' 2 ) 



A h\<s> 2 ) = Ah\<$> 2 ) 

= hA\<$> 2 ) 

= Ha 2 \§ 2 ) 

= a 2 H\§ 2 ) 



Where we have used the given that AH = HA. We may now write 

A\%)=a 2 \%) 



(4.5) 



(4.6) 
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Now consider again equation (4.1). If we take the adjoint of this equation we 
obtain 

($1^ = <$iK (4-7) 

Now use the fact that A = (we assumed A was Hermitian) and that the 
eigenvalues of a Hermitian operator are real. This yields 

($i|A= ($i|ai (4.8) 

Multiply on the right by |$' 2 ) 

($ 1 |i|$ , 2 ) = a 1 ($ 1 |$ , 2 ) (4.9) 
Now multiply equation (4.6) on the left by ($i| to obtain 

($ 1 |i|$ , 2 > = a 2 <$ 1 |$ , 2 ) (4.10) 

If we subtract equation (4.10) from equation (4.9) we arrive at 

{a 1 -a 2 )(<!> 1 \<!>' 2 )=0 (4.11) 

Since we assumed a\ ^ a 2 , then ($i|$ 2 ) = 0. Recalling the definition of |$ 2 ), we 
have 

($i|#|<I> 2 } = (4.12) 

which was to be proven. Thus if our desired wavefunction is an eigenfunction of 
some Hermitian operator that commutes with the Hamiltonian, our CI space need 
not include those iV-electron functions which are eigenfunctions of this operator 
with different eigenvalues. As an example, consider the spin angular momentum 
operator 5 2 . If we want to solve for a state of spin S, then we know 

5 2 |*) = 5(5+l)|*> (4.13) 

and any basis function of a different spin can be excluded from the CI. Slater 
determinants are generally not eigenfunctions of 5 2 . However, we can take linear 
combinations of Slater determinants so that we do have eigenfunctions of 5 , 
such functions are generally called configuration state functions, or CSF's. The 
advantage of using CSF's is that we can throw out all functions with the wrong 
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eigenvalue S — they contribute to another, noninteracting block of the Hamiltonian 
matrix. This reasoning also applies to symmetry operations of point groups, such 
that we can throw away any iV-electron basis functions (whether determinants or 
CSF's) which have the wrong irreducible representation. We can also restrict the 
basis functions according to their eigenvalues with respect to the operator S z . For 
a triplet state, we can perform the calculation using basis functions which have 
M s = —1,0, or 1. If we were to include basis functions of all these values of M s , 
we would obtain a triply-degenerate answer — as one should expect! 3 

4.2 Classification of Basis Functions by Excitation Level 

Now we will discuss the importance of various excitation classes to the CI wave- 
function. As noted in equation (2.9), the CI expansion is typically truncated 
according to excitation level; in the vast majority of CI studies, the expansion is 
truncated (for computational tractability) at doubly-excited configurations. Since 
the Hamiltonian contains only two-body terms, only singles and doubles can in- 
teract directly with the reference (for the sake of simplicity, we are assuming only 
a single reference for now). This is a direct result of Slater's Rules (cf. section 
2.4). The structure of the CI matrix with respect to excitation level is given be- 
low (adapted from Szabo and Ostlund [1], p. 235), where \S), \D), |T), and \Q) 
represent blocks of singly, doubly, triply, and quadruply excited determinants, re- 
spectively. The Hamiltonian matrix H is Hermitian; if only real orbitals are used, 
as is usually the case, then the Hamiltonian is also symmetric. Thus only the 
lower triangle of H is shown below. 



($o|#|$o) 




(S\H\S) 



H = 





(4.14) 



<T| 
(Q\ 












(T\H\D) (T\H\T) 
(Q\H\D) (Q\H\T) (Q\H\Q) 



3 However, state-of-the-art determinant-based CI algorithms often include computational simplifications if the 
M s = component is used [16]. 
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Note that the matrix elements (S\H\&o) are given as 0. This is due to Bril- 
louin's theorem, which is valid when our reference function |$q) is obtained by the 
Hartree-Fock method (Hartree-Fock guarantees that off-diagonal elements of the 
Fock matrix are zero, and it turns out that the matrix element between two Slater 
determinants which differ by one spin orbital is equal to an off-diagonal element 
of the Fock matrix). Furthermore, the blocks (X\H\Y) which are not necessarily 
zero may still be sparse; for example, the matrix element ($^1 |#|$*cfe/ )? which 
belongs to the block (D\H\Q), will be nonzero only if a and b are contained in 
the set {c, d, e, /} and if r and s are contained in the set {t, u, v, w}. 

Since only the doubles interact directly with the Hartree-Fock reference, we 
expect double excitations to make the largest contributions to the CI wavefunc- 
tion, after the reference state. Indeed, this is what is observed. Even though 
singles, triples, etc. do not interact directly with the reference, they can still be- 
come part of the CI wavefunction (i.e. have non-zero coefficients) because they 
mix with the doubles, directly or indirectly. Although singles are much less im- 
portant to the energy than doubles, they are generally included in CI treatments 
because of their relatively small number and because of their greater importance 
in describing one-electron properties (dipole moment, etc.) 

4.3 Energy Contributions of the Various Excitation Levels 

Table 2 demonstrates the importance of various excitation classes in obtaining CI 
energies. We see that singles and doubles account for 95% of the correlation en- 
ergy at the equilibrium geometries of the molecules listed. We see that quadruple 
excitations are more important than triples, at least as far as the energy is con- 
cerned. At stretched geometries, the CISD and CISDT methods become markedly 
poorer, yet the CISDTQ method still recovers a very high (and nearly constant) 
fraction of the correlation energy, suggesting that CISDTQ should give reliable 
results for energy differences across potential energy surfaces for molecules of this 
size. 
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Table 2: Percentage of correlation energy recovered by various CI excitation levels for some small 
molecules. 



Molecule 


Percent Corr. 


Energy 


CISD 


CISDT 


CISDTQ 


BH 


94.91 


n/a 


99.97 


H 2 0(R e ) 


94.70 


95.47 


99.82 


H 2 0(1.5 R^) 


89.39 


91.15 


99.48 


H 2 O(2.0 R e ) 


80.51 


83.96 


98.60 


NH 3 


94.44 


95.43 


99.84 


HF 


95.41 


96.49 


99.86 


H| 


96.36 


96.87 


99.96 



"Data from reference [6] except for H^" data 
from reference [17]. 



4.4 Size of the CI Space as a Function of Excitation Level 

We can also see from Table 3 that the number of A-electron basis functions in- 
creases dramatically with increasing excitation level. It should be pointed out that 
while the calculations on BH, HF, and Hy" used DZP basis sets, those on H2O and 
NH3 used only DZ basis sets. A DZP basis should be considred the minimum ad- 
equate basis for a truly meaningful benchmark study, and even then it is desirable 
to use a high-quality basis such as an Atomic Natural Orbital (ANO) set. While 
it is generally possible to perform CISD calculations on small molecules with a 
good one-electron basis, the CISDTQ method is limited to molecules containing 
very few heavy atoms, due to the extreme number of A-electron basis functions 
required. Full CI calculations are of course even more difficult to perform, so 
that despite their importance as benchmarks, few full CI energies using flexible 
one-electron basis sets have been obtained. 

The size of the full CI space in CSF's can be calculated (including spin sym- 
metry but ignoring spatial symmetry) by Weyl's dimension formula as applied to 
the Distinct row table (DRT). If N is the number of electrons, n is the number 
of orbitals, and S is the total spin, then the dimension of the CI space in CSF's 
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CSF 


's required 61 




Molecule 


CISD 


CISDT 


CISDTQ 


FCI 


BH 


568 


n/a 


28 698 


132 686 


H 2 


361 


3 203 


17 678 


256 473 


NH 3 


461 


4 029 


19 925 


137 321 


HF 


552 


6 712 


48 963 


944 348 


H| 


1 271 


24 468 


248 149 


2 923 933 



"Data from reference [6] except for Hy" data from 
reference [17]. 



is given by 

25 + 1 / n + 1 \( n + 1 \ 
D " ns = ^+T{n/2-S ) [n/2 + S+I j (4 ' 15) 

The dimension of the full CI space in determinants (again, ignoring spatial 
symmetry) is computed simply by 



JW, = [ N J[ Np ) (4.16) 
or, in a form closer to equation 4.15, 

DnNS = ( N/2 + S ) ( N/2- S ) (4A7) 

Table 4 shows the dimension of the full CI space (neglecting spatial symmetry) in 
determinants and in CSF's. Current full CI algorithms are limited to a few million 
determinants. Although there have been reports [18, 19, 20] of larger calculations 
(including a few billion determinants), the computational expense is (currently) 
too great for routine calculations of this size. 
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Number of electrons 




Orbitals 


6 


8 


10 


12 


10 


14.4 x 10 3 
(4.95 x 10 3 ) 


44.1 x 10 3 
(13.9 x 10 3 ) 


63.5 x 10 3 
(19.4 x 10 3 ) 


44.1 x 10 3 
(13.9 x 10 3 ) 


20 


1.30 x 10 6 
(379 x 10 3 ) 


23.5 x 10 6 
(5.80 x 10 6 ) 


240 x 10 6 
(52.6 x 10 6 ) 


1.50 x 10 9 
(300 x 10 6 ) 


30 


16.5 x 10 6 
(4.56 x 10 6 ) 


751 x 10 6 
(172 x 10 6 ) 


20.3 x 10 9 
(4.04 x 10 9 ) 


353 x 10 9 
(62.5 x 10 9 ) 



4.5 The Frozen Core Approximation 

It is quite common in applications of the CI method to invoke the frozen core ap- 
proximation, in which the lowest-lying molecular orbitals (occupied by the inner- 
shell electrons) are constrained to remain doubly-occupied in all configurations. 
The frozen core for atoms lithium to neon typically consists of the Is atomic or- 
bital, while that for atoms sodium to argon consists of the atomic orbitals Is, 2s, 
2p x , 2p y and 2p z . The frozen molecular orbitals are those which are primarily 
these inner-shell atomic orbitals (or linear combinations thereof). 

A justification for this approximation is that the inner-shell electrons of an 
atom are less sensitive to their environment than are the valence electrons. Thus 
the error introduced by freezing the core orbitals is nearly constant for molecules 
containing the same types of atoms. In fact, it is sometimes recommended that one 
employ the frozen core approximation as a general rule because most of the basis 
sets commonly used in quantum chemical calculations do not provide sufficient 
flexibility in the core region to accurately describe the correlation of the core 
electrons. 

Not only does the frozen core approximation reduce the number of configura- 
tions in the CI procedure, but it also reduces the computational effort required 
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to evaluate matrix elements between the configurations which remain. Assum- 
ing that all frozen core orbitals are doubly occupied and orthogonal to all other 
molecular orbitals, then it can be shown [21] that 

($/|# |$j> = ($i\H \$j) (4.18) 

where $/ and <lj are identical to $/ and $j, respectively, except that the core 
orbitals have been deleted from <lj and and H has been replaced by Ho defined 
by 

N-N c N-N c 1 

Ho = E c + £ E — (4-19) 

i=l i>j fij 

where N is the number of electrons and N c is the number of core electrons. E c 
is the so-called "frozen-core energy," which is the expectation value of the deter- 
minant formed from only the N c core electrons doubly occupying the n c = N c /2 
core orbitals 

n c n c 

E c = 2^hii + ^ {2(n|j'i) - (ij\ji)} (4.20) 

i ij 

Finally, h c (i) is the one-electron Hamiltonian operator for electron i in the average 
field produced by the N c core electrons, 

h c (i) = h(i) + Z fa® - kM (4.21) 

with Jj(i) and Kj(i) representing the standard Coulomb and exchange operators, 
respectively. Note that, although we have written the frozen core energy E c 
and frozen core operator h c in terms of molecular orbitals, it is not necessary to 
explicitly transform the one- and two-electron integrals involving core orbitals. 
Assuming real orbitals, we can define a frozen core density matrix [22] in atomic 
(or symmetry) orbitals as 

p;* = tc;ci (4.22) 

i 

where C* is the contribution of atomic orbital p to molecular orbital i. Now the 
frozen core operator in atomic orbitals becomes 

Ku = V + 2 £(HH^ C , - Upvl^Ka (4-23) 
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and the frozen core operator in molecular orbitals hfj can be obtained simply by 
transforming h c ^ y . Similarly the frozen core energy can be evaluated as 

Ec = E^(V + ^,) (4-24) 
= Tr{P c h) + Tr{P c h c ) 



An analogous approximation is the deleted virtual approximation, whereby a 
few of the highest-lying virtual (unoccupied) molecular orbitals are constrained 
to remain unoccupied in all configurations. Since these orbitals can never be 
occupied, they can be removed from the CI procedure entirely because no terms 
involving them contribute to the CI coefficients or energy. The rationalization for 
this procedure is that it is unlikely that electrons will choose to partially populate 
high-energy orbitals in their attempt to avoid other electrons. However, this 
conclusion is generally true only for very high-lying virtual orbitals (such as those 
formed by antisymmetric combinations of symmetry orbitals in the core region). 
For all other virtual orbitals, such simplistic reasoning is not sufficient. 

Davidson points out that those high energy SCF virtual orbitals which result 
from the antisymmetric combination of the two basis functions describing each 
atomic orbital in a double-C basis set (such as the 3p-like orbital formed from the 
minus combination of the larger and smaller 2p atomic orbitals on oxygen) often 
make the largest contribution to the correlation energy in M0ller-Plesset (MPn) 
wavefunctions [23]. This can be seen from the expression for the second-order 
correction to the energy in M0ller-Plesset perturbation theory, 

E(2) = E E \(ab\\r S )\^ {4 25) 

a>br>s Ca ~T~ Cfe € r 6 S 

where €i is the orbital energy (eigenvalue) of orbital %. Thus a virtual orbital r 
with a large energy e r will contribute to a large energy denominator in each term 
of equation (4.25) in which it appears. However, if orbital r lies close spatially 
to one of the orbitals a or 6 occupied in the reference, then this large overlap 
will contribute to a large two-electron integral (ab\\rs). This integral is squared in 



4 REDUCING THE CI SPACE 



32 



the numerator, leading to a large overall contribution to the second-order energy. 
For the antisymmetric combinations of the basis functions describing the core 
region, this large numerator is insufficient to overcome the even larger energy 
denominator; such virtual orbitals can generally be deleted with minimal loss in 
the correlation energy recovered. 

Although the analysis in the preceeding paragraph is based on perturbation 
theory, similar conclusions can be drawn for the CI method. This is most eas- 
ily verified by actual calculations, since analytical expressions for the energetic 
contribution of orbitals to the CI energy are not nearly as simple to obtain or 
interpret as is equation (4.25). 
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4.6 Truncated CI is not Size Extensive 

As we have previously pointed out, full CI — being the matrix formulation of the 
Schrodinger equation — is an exact theory for nonrelativistic electronic structure 
problems. If we truncate the CI (either in the one-electron or iV-electron space), 
we no longer have an exact theory. Of course either of these truncations will 
introduce an error in the wavefunction, which will cause errors in the energy 
and all other properties. One particularly unwelcome result of truncating the N- 
electron basis is that the CI energies obtained are no longer size extensive or size 
consistent. 

These two terms, size extensive and size consistent, are used somewhat loosely 
in the literature. Of the two, size extensivity is the most well-defined. A method 
is said to be size extensive if the energy calculated thereby scales linearly with 
the number of particles N. The word "extensive" is used in the same sense as in 
thermodynamics, when we refer to an extensive, rather than an intensive property. 
A method is called size consistent if it gives an energy Ea + Eb for two well 
separated subsystems A and B. While the definition of size extensivity applies 
at any geometry, the concept of size consistency applies only in the limiting case 
of infinite separation. In addition, size consistency usually also implies correct 
dissociation into fragments; this is the source of much of the confusion arising 
from this term. Thus restricted Hartree-Fock (RHF) is size extensive, but it 
is not necessarily size consistent, since it cannot properly describe dissociation 
into open-shell fragments. It can be shown that many-body perturbation theory 
(MBPT) and coupled-cluster (CC) methods are size extensive, but they will be 
size consistent only if they are based on reference wavefunction which dissociates 
properly. 

As previously stated, truncated CPs are neither size extensive nor size consis- 
tent. A simple (and often used!) example is sufficient to make the point. Consider 
two noninteracting hydrogen molecules. If the CISD method is used, then the en- 
ergy of the two molecules at large separation will not be the same as the sum 
of their energies when calculated separately. In order for this to be the case, 
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we would have to include quadruple excitations in the supermolecule calculation, 
since local double excitations could happen simultaneously on A and B. 

We would tend to think that size extensivity and size consistency are important, 
physical properties that all quantum mechanical models should have (indeed, full 
CI, an exact theory, has these properties), but perhaps they are not as essential 
as all that. Duch and Diercksen have claimed that "making size extensivity the 
most important requirement of quantum chemical methods, although it does not 
guarantee correct physical description, seems to be based not that much on physi- 
cal as on esthetical criteria" [24]. Indeed, they show that quantum mechanics is a 
"holistic" theory, not well-suited toward the description of separated subsystems: 

Hilbert space of antisymmetric, many particle functions, describing 
the total system, can not be decomposed into separate subspaces. Con- 
sider two systems, Sa and Sb, with Na and iVg electrons, respectively. 
Each system is described by its own function, antisymmetric in Na 
particles and ^/^ in Nb- Assuming that both functions are normalized to 
unity it is easy to show that the product function ^ ab = ^ b is always 
"far" from the antisymmetric function \I/ = A$> ab-, as measured by the 
overlap ab\^) or the norm of the difference 2 — \/2 < \\^ab — ^|| 2 < 2. 

Such arguments notwithstanding, it is clear that the fraction of the correlation 
energy recovered by a truncated CI will diminish as the size of the system increases, 
making it a progressively less accurate method. There have been many attempts 
to correct the CI energy to make it size extensive. The most widely-used (and 
simplest) of these methods is referred to as the Davidson correction [25], which is 



where E$d is the basis set correlation energy recovered by a CISD procedure. 
This correction approximately accounts for the effects of "unlinked quadruple" 
excitations (i.e. simultaneous pairs of double excitations). The multireference 
version [24] of this correction is 




(4.26) 




(4.27) 
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where Emrci is the multireference CI energy and Emr is the energy obtained from 
the set of references (MCSCF energy if the references are obtained as all references 
in an MCSCF procedure). We have simply replaced the CISD correlation energy 
in equation (4.26) with the analogous multireference correlation energy, and we 
have replaced Cq with the analogous sum of squares of all the reference coefficients. 

There are a number of other size extensivity corrections, and most of them 
do not take any significant amount of computation. Reference [24] provides a 
nice comparison of several of the more common CI size extensivity correction 
methods. We should also mention that Malrieu and co-workers have presented a 
self-consistent dressing of the Hamiltonian which gives size extensive results for 
selected CI procedures [26]. 
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5 Second Quantization 

Much of the literature in CI theory makes use of the notation of second-quantization 
Szabo and Ostlund [1] give a good introduction to second-quantized operators. 
Here we will only summarize the anticommutation relations between creation and 
annihilation operators, and then proceed to express the Hamiltonian in second 
quantized form for spatial orbitals, rather than for spin orbitals. Then we will use 
these results to derive the Hamiltonian in terms of the unitary group generators. 

The anticommutation relations for two annihilation operators is 



The anticommutation relation between a creation and an annihilation operator is 



Now we will find an expression for the Hamiltonian in terms of creation and 
annihilation operators over spatial orbitals. We begin with the second-quantized 
form of the one- and two-electron operators (see Szabo and Ostlund [1] p. 95) 




(5.2) 



(5.1) 



{ai,a]} = a l a] + a]ai = % 



(5.3) 



6i = Y,(i\h\j)4a } 



(5.4) 



= r Y.{ij\kl)a\a)aia k 



(5.5) 



ijkl 



where the sums run over all spin orbitals {xi}- Thus the Hamiltonian is 



H = i>Ja g [p|%] + - E alala s a q [pq\rs] 



2n \ 2n 



(5.6) 



pq Z pqrs 



From the previous equation we can see that the second-quantized form of the 
Hamiltonian is independent of the number of electrons in the system. Now inte- 
grate over spin, assuming that spatial orbitals are constrained to be identical for 
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a and (3 spins. A sum over all 2n spin orbitals can be split up into two sums, one 
over n orbitals with a spin, and one over n orbitals with (3 spin. Symbolically, 
this is 

2n n n 

E = E + E (5.7) 

a a a 

The one-electron part of the Hamiltonian becomes 

n 

Hone = EN%H<*°V + \p\ h \q] a pa a qP + [pI%K/3<V + [£|%K/3 a <Z/3 (5.8) 
pq 

After integrating over spin, this becomes 

+ + 

I J I I 1 / 1 U " -I- /T» 



#one = E(p|%){4a°V + a ll3 a <lp} ( 5 - 9 ) 



The two-electron term can be expanded similarly to give 
1 



^two — n E (P#| rS ){ a j)a a ra a sa a tfa + a p a a rP a sP a qa + a \j3 a ra a 3a^ql3 + a p/3 a J/3 a s/3 a g/?} 

(5.10) 

Now we make use of the anticommutation relation (5.1) and we swap the order of 
a s and a q , introducing a minus sign. This yields 

Htwo = "„ E (pq\ r s){cip a al a a qa a sa + a pa a r ^a qa a s p + a p ^al a a q pa sa + a p pa r pa q pa s p} 

- pqrs 

(5.11) 

Now we use the anticommutation relation between a creation and an annihilation 
operator, which is given by (5.3). This relation allows us to swap the a q and a\ 
in each term, to give 

H-two = ~ E (p^|^s) O'pa^qa^'ra^'sa 3qa,raQ>pa®'sa ^pa^q^rfi^sP 3qa,rf3Q- 'pa^s \/3 
- pqrs ' 

+ a) p pa q f3a\, a a sa — S q p^ ra a p ^a sa + a p pa q pa r pa s p — S q p jr pa p pa s p (5.12) 

Now we observe that 8 qa ,ra and 5 q p :r p can both be written 5 qr , and also that 5 qa>r p 
and 5 q f3, ra are both 0. This simplifies our equation to 

1 n 

Htwo = o E (pq\rs) 

l^pa^QCt^ra^sa ^pa^qa^rP^sP ^pP^qP^ra^sc^ + a pP a qP a rP a &P 
£ pqrs L 



5 SECOND QUANTIZATION 



38 



3qr a pa a sa $qr a p j3 a s(3 (5.13) 

Now we introduce the replacement (or shift) operator 

E i:j = a} a a ja + aipcijp (5.14) 

which Paldus has shown [5] to be isomorphic to the generators of the unitary 
group. This replacement operator is commonly referred to as a unitary group 
generator, but as Duch has pointed out [27], such usage is somewhat dubious in 
papers where no unitary group theory is employed. 

n 1 n / \ 

H = Y.(p\h\q)E pq + - E (PQ\rs) (E pq E rs - 5 qr E ps ) (5.15) 

P q Z pqrs ^ ' 

This is the Hamiltonian in terms of replacement operators. Contemporary papers 
on CI theory often express the Hamiltonian in the form of equation (5. 15). 4 



4 Reference [28] contains a minor mistake, giving (pq\rs) as ( P q\rs). 
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6 Determinant-Based CI 

6.1 Introduction to Determinant CI 

As we have already pointed out, the size of the CI space can be reduced signif- 
icantly by including only those iV-electron basis functions which have the same 
value of the quantum number S as the desired approximate wavefunction (cf. 
sections 4.1 and 4.4). Thus it would seem that one should always prefer CSF's 
to Slater determinants when performing a CI. However, certain computational 
advantages arise from using determinants, as we will discuss in the next few sec- 
tions. Many modern algorithms for performing large CPs (i.e. more than single 
and double excitations) use determinants as the iV-electron basis. 

Handy's 1980 paper "Multi-Root Configuration Interaction Calculations" [29] 
represented a major advance in determinant-based CI, even though the paper was 
more concerned with how integrals and CI coefficients are stored than with the 
computational advantages of determinants over CSF's. Handy used the Cooper- 
Nesbet method for performing the CI iteration; the CI coefficients are updated 
according to the formula 

fc - = (A) (61) 

where the r vector is defined as 

n = - E5 ij) c j (6-2) 

i 

Handy realized that, if determinants are used as basis functions, and particularly 
if these determinants are expressed as "alpha strings" and "beta strings," then Hij 
(and thus r) could be evaluated very efficiently. In order to understand Handy's 
reasoning, we must first define alpha and beta strings. 
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6.2 Alpha and Beta Strings 

Although Handy was the first to use alpha and beta strings, we will use the 
notation of Olsen et al. [16]. We define an alpha string as an ordered product of 
creation operators for spin orbitals with alpha spin. If I a contains a list {i, j, . . . n} 
of the n a occupied spin orbitals with alpha spin in determinant \I), then the alpha 
string a(I a ) is a\ a a^ a . . . o) na . A beta string is defined similarly. Thus we can 
rewrite Slater determinant \I) in terms of alpha and beta strings. 

\I) = \a(I a )P(I p ))=a(I a )P(I p )\) (6.3) 

For instance, suppose we have the Slater determinant \I) = |0i a 02a03a0i/302/3^4/?)- 
Then the alpha string a(I a ) is given by 

a (I a ) = a\ a a\ a a\ a (6.4) 

and the beta string is given by 

Note that the order of the creation operators matters; if we swap the order of two 
creation operators within the alpha string (or within the beta string), then we 
introduce a sign change (see equation 5.2). Also, acting the alpha string on the 
vacuum first, rather than the beta string, may introduce a minus sign, depending 
on the number of alpha and beta electrons. Although the order of the orbitals 
and whether the alpha or beta string acts first is of no real consequence, we must 
be sure to keep our use of alpha and beta strings consistent, or sign problems will 
result. In most of the literature, and in these notes, the beta string will be placed 
to the right of the alpha string in equations like (6.3). Further, within each string, 
orbitals are listed in strictly increasing order. 

Handy realized the following advantages to alpha and beta strings: 

1. Direct CI methods often require an index vector which points to a list of all 
allowed excitations from a given iV-electron basis function. Using alpha and 
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beta strings, the index vector need not be the length of the CI vector — its size 
is dictated by the number of alpha or beta strings, which is approximately 
the square root of the number of determinants. This results from the fact 
that (in determinant-based CI) electrons in alpha spin-orbitals can be excited 
only to other alpha spin-orbitals, and electrons in beta spin-orbitals can be 
excited only to other beta spin-orbitals. 5 

2. To form r(I a ,Ip) in equation (6.2), all functions \a(J a )(3(Jp)} which have 
non-zero matrix elements with \a(I a )(3(Ip)) are generated, one at a time, 
with the appropriate integral being looked up and multiplied by the appro- 
priate CI coefficient. No time is wasted considering determinants which are 
noninteracting, and the coefficients of the integrals are simply ±1. 

3. Efficiency is increased by realizing that all integrals which enter the expression 
(a(I a )(3(Ip)\H\a(J a )(3(Ip)) (equation 6.2), where a(J a ) differs from a(I a ) by 
two orbitals, are independent of (3(1 p) ■ 

We will make these points more clear in our discussion of the RAS CI, which 
is a direct extension of Handy's observations concerning alpha and beta strings. 
However, at this point we will proceed to discuss the graphical representation of 
alpha and beta strings. 

6.2.1 Graphical Representation of Alpha and Beta Strings 

Any CI program requires some method for assigning a numerical code to each 
determinant (or configuration). Determinant-based CPs also require an address- 
ing scheme for each alpha and beta string. Such addressing can be facilitated by 
using graphical representations. This section discusses the construction of simple 
graphs representing alpha and beta strings, and how the numerical code or ad- 
dress of each string can be calculated from these graphs. Our approach will be 
based on the work of Duch, who has described [27] the graphical representation of 
CI spaces in considerable detail. First, we consider the simple two-slope directed 

5 In the determinant CI expansion, we restrict all determinants to a single value of M s . 
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graphs ("digraphs") Qiip, : n a /p) which represent alpha or beta strings without 
consideration of point-group symmetry. Figure 1 presents such a digraph repre- 
senting all strings with five electrons in seven orbitals (as might be appropriate 
for the alpha or beta strings of H2O in a minimal basis set, neglecting point-group 
symmetry). Each string is represented by a "walk" on the graph, from the head 
(at e = o = 0) to the tail (at e = n a /p, o = n); the graph should always be 
traversed in this direction. Moving straight down from vertex (e, o) to vertex 
(e, o + l) indicates that orbital o + 1 is unoccupied in the current string, while 
moving down diagonally from vertex (e,o) to vertex (e + l,o + 1) indicates that 
orbital o + 1 is occupied. Each vertex on the graph is assigned a weight x(e,o), 
and each arc connecting two vertices is assigned an arc weight Y(e,o) for the arc 
leaving vertex (e, o). Since, in general, two different arcs can leave a given vertex, 
we write Yo(e, o) for the arc originating from vertex (e, o) which leaves orbital o + l 
unoccupied, and Yi(e,o) for the arc which occupies orbital 0+ l. 6 The index or 
address of a string or walk is obtained by adding weights for each arc contained 
in the walk, i.e. 

n 

I a {L a ) = X{L a ) + Y.y L X^) (6-6) 

where Li is the occupation (0 or 1) of the zth arc, and [e^i) are the coordinates 
of the vertices crossed by L a . The term X(L a ) gives the offset of a given graph, 
if more than one graph is employed. The index for a determinant is given by 
I(L a ,L,P) = I a (L a )S a + Ip(L^), where S a is the number of alpha strings. This 
leads us to write the CI coefficients as a rectangular matrix C(I a , Ip). Restrictions 
on the CI space mean that only certain subblocks of the C matrix are allowed to 
be nonzero. 

There are several different methods for assigning the arc weights by which we 
calculate the index of a string according to equation (6.6). Under the lexical 
ordering scheme, the tail (n a , n) of an alpha string graph is assigned a weight 
x = 1. Other vertex weights are computed according to the recursive formula 

x(e, o) = x(e + 1, o + 1) + x(e, o + 1) (6.7) 

6 This differs somewhat from Duch [27], who sometimes uses Y(e,o) to denote the arc entering vertex (e, o) in 
reverse-lexical addressing. 
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Using lexical ordering, typically all arc weights Yo(e, o) are set equal to zero, and 
the arc weights Yi(e,o) are determined according to 

Yi(e, o) = x(e + 1, o + 1) + x{e + 1, o) H h x(e + 1, e + 1) (6.8) 

Figure la features vertex and arc weights computed in this manner. A result of 
the lexical ordering scheme is that paths with a fixed upper part and an arbitrary 
lower part have contiguous indices. The particular choice of Y values above is 
appropriate if the rightmost path is to have an index of zero. The same effects 
can be achieved using 

Y 1 {e : o) = (6.9) 
y (c,o) = x(e + l,o+l) (6.10) 

as illustrated in Figure lb. Any walk has the same index in Figures la and lb. 
For instance, the walk a2 a a\ a a\ a a\ a a\ a has an index of 5 + 4 + 3 + 2 + 2 = 16 
(equation 6.6) from Figure la, and an index of 15 + 1 = 16 from Figure lb. 

In the so-called "reverse-lexical" ordering scheme, all upper paths for a fixed 
lower path have contiguous indices. Vertex weights are now determined as 

x(e, o) = x(e, o — 1) + x(e — 1, o — 1) (6-11) 

where the overbar indicates reversed-lexical ordering. Figure 2a depicts a reversed- 
lexical graph with all non-occupied orbital arcs set to zero. The occupied orbital 
arcs are computed as 

Yi(e,o) = x(e,o) (6.12) 

Figure 2b is the same except that now all occupied arcs have weights of zero. The 
non-occupied arc weights are 

?o(e, 6) = x(e, 6) + x(e + 1, o + 1) H h x{N -l,o + AT-e-l) (6.13) 

Note that string indices for reverse-lexical ordering are not necessarily the same as 
indices for lexical ordering. For the string al a al a a\ a a\ a a\ a considered previously, 
the index is calculated as 1 + 1 + 1 + 1+ 6 = 10 from Figure 2a, or as 5 + 5 = 10 
from Figure 2b. 
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Figure 1: Alpha string graph for n a = 5, n = 7. Vertex weights are determined according to lexical 
ordering, and arc weights are given so that the rightmost path has index zero, (a) All unoccupied 
arc weights Y (e,o) are zero, (b) All occupied arc weights Yi(e, o) are zero. 



e e 
012345 012345 

i i i i i i I i i i i i r 




(a) 



(b) 



6 DETERMINANT CI 



45 



Figure 2: Alpha string graph for n a = 5,n = 7. Vertex weights are determined according to 
reverse-lexical ordering, and arc weights are given so that the rightmost path has index zero, (a) 
All unoccupied arc weights Y (e,o) are zero, (b) Occupied arc weights Yi(e,o) are set to zero. 



e e 
012345 012345 

i i i i i i I i i i i i r 
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The arc weights given in Figures 1 and 2 cause the rightmost path to have an 
index I(R m ) = 0. If we change the arc weights so that the leftmost path has index 
I(L m ) = 0, we obtain four more addressing schemes. The two simplest schemes 
for I(L m ) = are 

y (e,o) = Y 1 (e,o)=x(e,o + l) (6.14) 
?i(e,o) = Y (e,o) = x(e-l,o) (6.15) 

where the overbars indicate that reversed-lexical vertex weights have been used. 
Alpha strings for 5 electrons in 7 orbitals employing these addressing schemes are 
depicted in Figure 3. 

If we add another coordinate T to each vertex, we can extend these simple 
digraphs to include point-group symmetry. 

6.3 Restricted Active Space CI 

The Restricted Active Space (RAS) CI method was introduced by Olsen, Roos, 
J0rgensen, and Aa. Jensen [16] in 1988. The RAS method calls for the partitioning 
of the one-electron basis into four subsets. The first subset consists of the core 
orbitals, which are constrained to remain doubly-occupied. The remaining three 
subsets are labelled I, II, and III, and the CI space is limited by requiring a 
minimum of p electrons in RAS I and a maximum of q electrons in RAS III. 
There are no restrictions on the number of electrons in RAS II, and thus it is 
analogous to the complete active space (CAS). The full CI can be obtained as the 
maximum limit of the RAS space. Interestingly, although the main focus of the 
paper is on the utility of the RAS method in limiting the size of CI calculations, 
the maximum impact of this paper has been on the development of determinant- 
based full CI algorithms [19, 20]. 

The RAS CI method relies on Handy's separation of determinants into alpha 
and beta strings (see section 6.2). As in other determinant-based CI methods, 
the basis determinants are restricted to those having a given value of M s . Since 
the number of electrons N is also fixed, this means that the alpha and beta 
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Figure 3: Alpha string graph for n a = 5,n = 7, with arc weights determined so that the leftmost 
path has index zero, (a) Vertex weights for lexical ordering, and arc weights according to Y (e, o) = 
0, Yi(e, 6) = x(e, o + 1). (b) Vertex weights according to reverse-lexical ordering, and arc weights 
according to Yi(e, 6) = 0, Y (e, o) = x(e — 1,6). 
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strings always have constant lengths of n a and np, respectively. For a full CI, 
one forms all possible alpha and beta strings for a given n a and np, and the basis 
determinants are all possible combinations of these alpha and beta strings. In a 
RAS CI, the CI space is restricted in two ways: first, not all alpha and beta strings 
are allowed, and secondly, not all combinations of alpha and beta strings to form 
determinants are accepted. This is best seen from an example: consider the case 
of 6 orbitals, with n a = np = 3. If orbitals 4, 5, and 6 constitute RAS III, with 
a maximum of 2 electrons allowed, then clearly alpha strings such as a\ a a\ a a\ a 
are not allowed. Similarly, even though a\ a a\ a a\ a and a\pa\pa\p are allowed alpha 
and beta strings, these strings cannot be combined with each other because the 
resulting determinant would represent a quadruple excition into RAS III. 

If we employ a graphical description of the alpha and beta strings as described 
in section 6.2.1, in general we require one graph for alpha strings and one graph for 
beta strings. However, in the case that M s = 0, only one graph is needed because 
the alpha string and beta string graphs are identical. As previously mentioned, 
for a RAS space not all alpha and beta strings can be freely combined. While it 
would be possible to create a table listing all allowed combinations of alpha and 
beta strings, there is a more efficient way around this difficulty Instead of using a 
single graph to represent all alpha (or beta) strings, instead we use several graphs. 
For example, we might use one graph for all strings with no electrons in RAS III, 
one graph for all strings with one electron in RAS III, and one graph for all strings 
with two electrons in RAS III. In this case, the restrictions on combinations of 
strings become restrictions on combinations of graphs — a more efficient treatment 
computationally. Figure 4 displays string graphs for n a = np = 3 and n = 6 for at 
most 2 electrons in RAS III, orbitals 4-6. Graph (a) represents all walks with two 
electrons in RAS III; graph (b) gives all walks with one electron in RAS III; and 
graph (c) gives the one walk with no electrons in RAS III. If only two electrons 
are allowed in RAS III, it is clear that alpha strings of graph a can be combined 
only with the beta string from graph (c), and alpha strings of graph (b) can be 
combined with beta strings of graphs (b) and (c). An alpha string from graph (c) 
can be combined with beta strings from graphs (a), (b), and (c). 
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Figure 4: String graphs for n a /p = 3,n = 6 with at most two electrons in RAS III of orbitals 4-6. 
Vertex weights and arc weights are given for lexical ordering. 
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6.3.1 Olsen's Full CI a Equations 

We now turn our attention to the general formulation of the full CI problem in 
terms of alpha and beta strings. Later on, we will consider how to modify our 
results for RAS CPs. We begin by describing Olsen's expressions for a, which is 
the action of the Hamiltonian H on the CI vector (or matrix) C(I a ,Ip). First we 
express H in second-quantized form (see section 5). 

H = Y,hkiE k i + - £fa|M) (EijEu - SjkEu) (6.16) 
hi 1 ijki v 

where E^i is the shift operator 

Em = E%i + E kl 

Again, a is given by 

<r(I a ,l0)= E (P(Jp)a(J a )\H\a(I a )(3(If3))C(J a , Jp) (6.18) 

As we will proceed to show, a can be split up into three terms: one involving only 
beta components of the linear group generators (o"i), one involving only alpha 
components of the generators (02), and one involving mixtures of the two (03). 
Inserting equation (6.16) into equation (6.18) yields 



(6.17) 
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= E (f3(Jp)a(J a )\Y,h k iE M (6.19) 

J ex. i J (3 kl 

+ \ t(ij\M) (EijEki ~ SjkEu) \a(I a )P(I ))C(J a , Jp) 
1 ijki v 7 

Now expanding the shift operators according to equation (6.17), we write a as 
a sum of three terms 

a(I a , Ip) = ai(I a , Ip) + a 2 (I a , Ip) + a 3 (I a , Ip) (6.20) 

where 

vi(IaJp) = E(W)«WIE^£ (6-21) 

J oiiJ kl 

+ \ tmO [HMi - \a(I a )(3(Ip))C(J a , Jp) 



2 ijki 



and 



n 



J J ft kl 

+ \ tmO {E-M - 8 jk E%) \a(I a )P(Ip))C(J a , Jp) 
z ijki v y 



and 

1 71 

°*{i a ,b) = E (W)«(4)kE(y'l^) + 4^S) MW(//0>c(j a , 



J a ,Jp * ijki 



(6.23) 



Obviously, these expressions can be simplified further. First, we observe that 
the expression for o\ contains no a operators. Therefore, (Ti(I a ,Ip) = unless 
J a = I a - Using this fact, and integrating out the a part, we obtain 



vi(iaJp) = UP(Jp)\Eh k iE p kl + - E«l*0 \%Ki - W) Wp))c(i a , jp) 



J f3 kl z ijki 



(6.24) 
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Taking out the Kronecker delta term and rearranging, we have 

Ti 

°i{IaM = Y.Y.hkmJp)\E P kl Wp))C{Ia,Jp) (6.25) 

J/3 kl 
1 n 

+ ~Y.Y,{m)Wp)\^KiWp))c{i«,jp) 

Z J/3 ijkl 

Now the second term can be combined with the first to give equation (9a) of 
reference [16]. 



Jfj kl 
1 

2 



1 n 

hki - xE(MiO 

Z 3 



C(I a ,J ) (6.26) 



Jp ijkl 

Similarly, o<i can be simplified to 

1 n 

C(J a ,I fi ) (6.27) 



°*{i a M = Y,Y,HJ*)\E a kl Wa)) 



j n 

hki - 9 E(MiO 



l n 

EE 



For efficient implementation, it is convenient to precompute the quantities 

^ n 

Ki = h kl - -Ete'b'O (6-28) 
Finally, we simplify 03. It may be rewritten as 



J Q ,J/3 ijfc/ 



+ o E E(/?(^)«(Ja)|^g^J«(/ a )/3(^))^1/c/)C(J a ,^) 



J a ,Jp ijkl 



(6.29) 
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Since we sum over all ijkl, we can permute i and j with k and /. We can also swap 
Efj and as can easily be verified. This yields equation (9c) from reference [16]. 

°*{IaM = E f:</3(^)|^g|/?(/ /3 )><o(^)|^ z |c.(/ a ))(zj>0^(^,^) 

JayJp ijkl 

(6.30) 

Thus we have written the action of the Hamiltonian on the current CI vector in 
terms of alpha and beta strings and alpha and beta shift operators. The product a 
is written as a sum of three terms: the first (<7i) involves only beta shift operators, 
the second (02) involves only alpha shift operators, and the third (03) involves 
both alpha and beta shift operators. Note that, except for the factor C(I a , Jp), 
o\ is independent of I a so that the algorithm for computing o\ is vectorizable 
(see below). Analogous results hold for o<i . This situation does not obtain in the 
computation of 03, however, which is the rate-limiting step. 

6.4 Full CI Algorithm 

It is easy to vectorize the formation of o\ since each element <Ji(I a: Ip) can be 
written as independent of I a apart from multiplication by a factor C(I a , Jp). The 
vectorized algorithm for the evaluation of 01, adapted from reference [16], appears 
in Figure 5. An analogous algorithm can be used to obtain 02. However, we can 
also obtain 02 for M s = states by 

a 2 (I a: I p ) = (-l) s a 1 (I p J a ) (6.31) 

as proven by Olsen and co-workers [16]. 

A fairly straightforward algorithm for evaluating 03 is presented in Figure 6. 
The vectorized algorithm for the evaluation of 03 is presented in Figure 7. This 
algorithm makes use of a gather and scatter operation to avoid indirect addressing. 

For M s = 0, an improvement to the 03 algorithm can be made by utilizing an 
identity similar to equation (6.31). The ijklth component of 03 is related to the 



DETERMINANT CI 



Figure 5: Vectorized Algorithm for U\. 

loop over Ip 

Set array F(Jp) = 

Loop over excitations E^ from \(3{Ip)) 

\P(K P )) = sgn(M)^|/3(7 /3 )) 

F(Kp)+ = sgn(kl)h' kl 

Loop over excitations E^ from \(3{Kp)) 

\/3{Jp))=s e a(ij)Et\l3(Kf l )) 
F(Jp)+ = (l/2)sgn(kl)sgn(ij)(ij\kl) 
end loop over E^ 
end loop over E kl 

ai(I a , Ip) = Y,j f3 F(J p )C{I a , Jp); vect'd over I a 
end loop over Ip 



Figure 6: Simple Algorithm for a 3 . 

loop over I a 

loop over \a(J a )) = sgn(kl)E^\a(I a )) 
loop over Ip 

loop over \(3(Jp)) = sgn{ij)E^{Ip)) 

a 3 (I a ,Ip) + = sgn(ij)sgn(kl)(ij\kl)C(J a , Jp) 
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Figure 7: Vectorized Algorithm for 0-3. 

loop over kl 

set up L(I), R(I), and sgn(J), defined by 
\a[L{I)]) = E%\a[R(l)])8ga(T) 
C'(I, Jp) = C[L(I), J /3 ]sgn(J); vectorized gather 
loop over Ip 

loop over excitations E^ from \/3(Ip)) 
\(3{Jp)) = sgn{ij)E^\(3{Ip)) 
F(Jp)+ = Sga (ij)(ij\kl) 
end loop over E^ 

V(I) = Ej, F(Jp)C'(I, Jp); vect'd over / 
a 3 [R(I), Ip]+ = V(J); vectorized scatter 
end loop over Ip 
end loop over kl 



klijih component via 

af\l a Jp) = (-1)^(^,4) (6-32) 
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